%Question 2 - Portfolio Optimization with a Twist
%Lu Xin lx108

function [ ] = Question2( )

clear all; clc; close all; warning('off');
% ---------------------------- SETUP ---------------------------- %
% Monthly Returns:
mu = [0.006 0.01 0.014 0.018 0.022];
% Correlation Matrix:
sig = [0.085 0.08 0.095 0.09 0.1];

corr = ones(5)*0.3 + eye(5)*0.7;
covMatrix = corr2cov(sig, corr);

% Short-Selling not allowed
% Returns corresponding to different months are mutually independent
% ---------------------------------------------------------------- %
% Call part 2a
% Get min variance
e = ones (1,5);
b = 1;
f = zeros (1,5);
[minVarWeights, pstdev] = markowitz(covMatrix, f, e, b);
rmin = mu * minVarWeights;

% Get max return
f = -mu;
H = zeros(5);
[maxRetWeights, ret] = markowitz(H, f, e, b);
rmax = -ret;

% ---------------------------------------------------------------- %
% Call part 2b
[trueFrontier, trueWeights] = optimization(covMatrix, mu, rmin, rmax);

%subplot(2,2,1);
%plot(trueFrontier(:,1),trueFrontier(:,2))
xlabel('Std Dev');
ylabel('Expected Return');
title('True Efficient Frontier');
% ---------------------------------------------------------------- %
% Call part 2c
[estFrontier, estWeights] = randomReturns(mu, covMatrix, 24);
%subplot(2,2,2);
%plot(estFrontier(:,1),estFrontier(:,2));
xlabel('Std Dev');
ylabel('Expected Return');
title('Estimated Efficient Frontier');

% ---------------------------------------------------------------- %
% Part 2d
[actualStdDev, actualReturns] = actualFrontier(estWeights, mu, covMatrix);
%subplot(2,2,3);
%plot(actualStdDev, actualReturns);
xlabel('Std Dev');
ylabel('Expected Return');
title('Actual Efficient Frontier');
% ---------------------------------------------------------------- %
% Part 2e
[aveExpectedFrontier, aveActualFrontier] = averagedFrontier(trueFrontier, ...
    mu, covMatrix, 1800, 10000);
% ---------------------------------------------------------------- %
end

% ---------------------------------------------------------------- %
% 2.a: Horizon of 1, markowitz model
function [weights, result] = markowitz(H, f, A, b)

lb = zeros(1, 5);
opts = optimset('Algorithm','active-set','Display','off');

[weights,result,exitflag,output,lambda] = ...
    quadprog(H,f,[],[],A,b,lb,[],[],opts);
end
% ---------------------------------------------------------------- %

% 2.b Plot efficient frontier
function [effFrontier, effWeights] = optimization(H, mu, rmin, rmax)
effFrontier = ones(10,2);
effWeights = ones (10,5);
f = zeros (1,5);
A = [mu; ones(1,5)];

for i=0:9
    rp = rmin + i * ((rmax - rmin)/9);
    b = [rp 1];
    [weights, result]  = markowitz(H, f, A, b);
    result  = sqrt(result * 2);
    effFrontier(i+1, :) = [result rp];
    effWeights(i+1, :) = weights;
end
end
% ---------------------------------------------------------------- %
% 2.c Random 
function [estiFrontier, estWeights] = randomReturns(mu, covMatrix, months)
randomReturns = mvnrnd(mu, covMatrix, months);
estMu = mean(randomReturns);
% sig = var(randomReturns);
covMatrix = cov(randomReturns);

% Get min variance portfolio
e = ones (1,5);
b = 1;

f = zeros (1,5);
[weights, pvar] = markowitz(covMatrix, f, e, b);
rmin = estMu * weights;

% Get max return portfolio
f = -estMu;
H = zeros(5);
[weights, ret] = markowitz(H, f, e, b);
rmax = -ret;

% Get Frontier and Weights
[estiFrontier, estWeights] = optimization(covMatrix, estMu, rmin, rmax);
end
% ---------------------------------------------------------------- %
% 2d - Actual
function [actualStdDev, actualReturns] = ...
    actualFrontier(estWeights, mu, covMatrix)
actualReturns = estWeights * mu';
actualStdDev = sqrt(diag(estWeights * covMatrix * estWeights'));
end
% ---------------------------------------------------------------- %
function [aveExpectedFrontier, aveActualFrontier] = averagedFrontier(trueFrontier, mu, covMatrix, months, runs)
sumExpectedFrontiers = zeros (10,2);
sumActualFrontiers = zeros (10,2);
% averagedFrontier
for j = 1:runs
    j
    [estFrontier, estWeights] = randomReturns(mu, covMatrix, months);
    sumExpectedFrontiers = sumExpectedFrontiers + estFrontier;
    
    [actualStdDev, actualReturns] = actualFrontier(estWeights, mu, covMatrix);
    sumActualFrontiers = sumActualFrontiers + [actualStdDev, actualReturns];
end

aveExpectedFrontier = sumExpectedFrontiers/runs;
aveActualFrontier = sumActualFrontiers/runs;
aveExpectedFrontier
aveActualFrontier
trueFrontier
%subplot (2,2,4);
plot(aveActualFrontier(:,1), aveActualFrontier(:,2),...
    aveExpectedFrontier(:,1),aveExpectedFrontier(:,2),...
    trueFrontier(:,1), trueFrontier(:,2));
xlabel('Std Dev');
ylabel('Expected Return of Portfolio');
legend('Actual Frontier','Expected Frontier','True Frontier');
title('10000 Runs, 1800 Months - Expected vs Actual vs True Efficient Frontier');
end

